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Abstract 

Shape  Memory  Alloys  (SMAs)  have  recently  been  considered  for  dynamic  loading 
applications  for  energy  absorbing  and  vibration  damping  devices.  An  SMA  body 
subjected  to  external  dynamic  loading  will  experience  large  inelastic  deformations 
that  will  propagate  through  the  body  as  phase  transformation  and/or  detwinning 
shock  waves.  The  wave  propagation  problem  in  a  cylindrical  polycrystalline  SMA 
rod  induced  by  an  impact  loading  is  considered  in  this  paper.  Numerical  solutions 
for  various  boundary  conditions  are  presented  for  stress  induced  martensite  and 
detwinning  of  martensite.  The  numerical  simulations  utilize  an  adaptive  Finite  El¬ 
ement  Method  (FEM)  based  on  the  Zienkiewicz-Zhu  (ZZ)  error  estimator.  Selected 
results  are  compared  to  known  analytical  solutions  to  verify  the  adaptive  FEM  ap¬ 
proach.  The  energy  dissipation  in  an  SMA  rod  is  evaluated  for  a  square  pulse  stress 
input  applied  at  various  temperatures  involving  both  stress  induced  martensite  and 
detwinning  of  martensite.  The  dynamic  response  of  a  NiTi  SMA  rod  is  also  studied 
experimentally  in  a  split  Hopkinson  bar  apparatus  under  detwinning  conditions. 
Strain  history  records  obtained  by  strain  gauges  placed  at  different  locations  along 
the  SMA  rod  are  compared  with  numerical  simulations  for  a  square  pulse  stress 
input.  The  quasi-static  and  dynamic  stress-strain  hysteretic  response  of  the  SMA, 
both  due  to  detwinning,  are  found  to  be  nearly  identical.  The  quasi-static  tests  are 
used  to  calibrate  the  rate  independent  constitutive  model  used  for  the  numerical 
simulations,  which  are  found  to  match  the  experimental  observations  reasonably 
well. 
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1  Introduction 


There  are  many  areas  of  applications  which  can  successfully  utilize  the  unique 
properties  of  SMAs.  The  engineering  research  presented  in  this  paper  relates 
directly  to  the  design  of  SMA  components  capable  of  absorbing  dynamic  loads. 
Such  components  can  be  integrated  into  critical  parts  of  structures  that  may 
need  protection  from  impact  loads.  Examples  include  joints  that  connect  the 
hull  of  an  underwater  vehicle  with  its  internal  structure,  tank  armor  or  blast 
resistant  cargo  containers.  Another  promising  field  of  application  includes  var¬ 
ious  active  or  passive  vibration  damping  devices.  Many  different  SMA  devices 
have  been  proposed  among  which  nonlinear  hysteretic  SMA  springs  (Yiu  and 
Regelbrugge,  1995;  Graesser,  1995),  wires  (Thomson  et  al.,  1995;  Fosdick  and 
Ketema,  1998)  or  rods  (Feng  and  Li,  1996).  In  a  recent  paper  (Lagoudas  et  ah, 
2001)  the  authors  investigate  numerically  the  vibration  damping  capabilities 
of  SMAs. 

Shape  Memory  Alloys  are  a  class  of  materials  that  change  their  internal  struc¬ 
ture  due  to  changes  in  temperature  and/or  externally  applied  loads.  At  high 
temperatures  the  crystal  lattice  is  in  the  high  symmetry  austenite  phase  (A). 
At  low  temperatures  the  material  exists  in  a  low  symmetry  martensite  phase 
(M).  The  austenite  to  martensite  phase  transition  is  diffusionless  and  is  char¬ 
acterized  by  shear  deformations  of  entire  regions  inside  the  material  (Way- 
man,  1983).  What  makes  SMA  materials  remarkably  different  from  ordinary 
metals  is  the  shape  memory  e^ect  and  the  effect  of  pseudoelasticity  which 
are  associated  with  the  specific  v/ay  the  phase  transition  occurs  (Funakubo, 
1987).  The  shape  memory  effect  allows  material  which  has  been  deformed 
while  in  the  martensitic  phase  to  recover  its  shape  upon  heating.  The  mech¬ 
anism  behind  this  behavior  is  the  ability  of  SMAs  to  allow  detwinning  of  the 
self-accomodated  martensitic  variants.  The  pseudelasticity  in  SMAs  is  their 
ability  to  support  large  inelastic  strains  recoverable  upon  unloading  due  to 
the  reverse  phase  transformation  from  martensite  into  austenite.  The  primary 
way  in  which  such  strains  are  introduced  in  the  material  is  the  stress  induced 
phase  transformation  from  austenite  into  martensite.  The  pseudoelastic  re¬ 
sponse  provides  both  energy  dissipation  capabilities  and  shape  recovery  dur¬ 
ing  the  thermomechanical  loading  path.  Utilizing  the  shape  memory  effect  also 
leads  to  dissipation  of  mechanical  energy  but  the  SMA  has  to  be  heated  after 
the  loading  is  applied  to  recover  its  shape. 

Several  constitutive  models  have  been  developed  in  recent  years  to  model  the 
shape  memory  effect  and  pseudoelasticity  of  polycrystalline  SMAs.  Among  the 
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most  widely  accepted  rate  independent  models  are  the  exponential  (Tanaka, 
1986),  cosine  (Liang  and  Rogers,  1990)  and  polynomial  (Boyd  and  Lagoudas, 
1994, 1996).  Any  of  these  models  can  be  unified  using  a  thermodynamic  frame¬ 
work  (Lagoudas  et  ah,  1996)  based  on  the  selection  of  appropriate  thermo¬ 
dynamic  potentials.  This  unified  constitutive  model  is  extended  bj'  Bo  and 
Lagoudas  (1999a,b,c,d)  to  incorporate  transformation  induced  plastic  defor¬ 
mations  and  to  account  for  the  evolution  of  the  material  behavior  during 
cyclic  loading.  In  the  models  proposed  by  Brinson  (1993);  Brinson  and  Lam- 
mering  (1993);  Lagoudas  and  Shu  (1999)  the  martensitic  volume  fraction  is 
subdivided  in  two  parts  to  account  for  thermally  induced  self-accommodated 
martensite  and  stress  induced  detwinned  martensite.  A  different  approach  is 
taken  by  Abeyaratne  et  al.  (1993, 1994);  Abeyaratne  and  Knowles  (1993)  who 
consider  a  rate  dependent  constitutive  model  that  allows  for  softening  duiing 
phase  transformation.  Other  authors  such  as  Patoor  et  al.  (1996),  Siderey  et  al. 
(1999);  Sun  and  Hwang  (1993a,b)  use  micromechanical  techniques  to  average 
the  response  of  the  parent  austenitic  phase  and  the  different  martensitic  vari¬ 
ants  to  obtain  a  model  for  the  macroscopic  behavior  of  polycrystalline  SMAs. 
For  further  details  on  SMA  models  the  reader  is  referred  to  the  work  by  Qidwai 
and  Lagoudas  (2000b)  as  well  as  the  review  paper  by  Birman  (1996).  In  the 
current  work  the  unified  approach  (Lagoudas  et  ah,  1996)  is  chosen  over  the 
more  complex  micromechanical  models,  assuming  rate  independence  in  the 
constitutive  thermomechanical  response  of  SMAs.  As  it  will  be  shown  later 
such  an  assumption  is  confirmed  experimentally  for  the  case  of  detw  inning  of 
martensite. 

In  a  recent  paper  (Chen  and  Lagoudas,  2000)  the  rate  independent  model  for 
polycrystalline  SMAs  (Lagoudas  et  al.,  1996)  is  employed  to  obtain  solutions 
to  the  coupled  thermomechanical  problem  for  SMA  materials.  The  authors 
take  into  account  the  latent  heat  generation  and  assuming  adiabatic  condi¬ 
tions  they  solve  the  problem  by  the  method  of  characteristics  together  with 
jump  conditions  that  yield  unique  solutions.  A  similar  study  (Bekker  et  al., 
2002),  but  for  different  constitutive  models  has  been  carried  out  for  both 
isothermal  and  adiabatic  conditions.  In  a  different  setting  Oberaigner  et  al. 
(1996)  investigates  numerically  the  coupled  problem  of  wave  propagation  and 
heat  transfer  in  an  SMA  rod.  The  authors  focus  on  stress  pulses  of  low  mag¬ 
nitude  that  cause  only  elastic  deformations.  The  temperature  at  one  end  of 
the  SMA  rod  is  chosen  as  a  function  of  time  in  a  such  a  way  as  to  utilize 
the  phase  change  due  to  the  shape  memory  effect  in  order  to  maximize  the 
damping  characteristics  of  the  rod. 

The  dynamics  of  phase  transformation  in  piecewise  linear  elastic  materials 
with  non-monotone  hysteresis  is  also  studied  by  Abeyaratne  and  Knowles 
(1991).  A  unique  solution  is  obtained  with  the  use  of  a  kinetic  relation  con¬ 
trolling  the  rate  of  the  phase  transformation  together  with  a  nucleation  con¬ 
dition  for  the  initiation  of  the  transformation.  In  later  work  the  same  authors 
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extend  the  analysis  to  account  for  thermal  effects  (Abeyaratne  and  Knowles, 
1994a, b).  In  a  general  setting  Pence  (1986)  considers  wave  propagating  in  a 
nonlinear  elastic  bar  with  a  non-monotonic  stress-strain  relationship  subjected 
to  a  monotonically  increasing  load.  It  is  found  that  for  sufficiently  high  loads 
a  strain  discontinuity  associated  with  phase  transformation  is  being  created. 

There  has  been  a  limited  amount  of  experimental  work  done  on  characterizing 
the  dynamic  response  of  SMAs.  An  experimental  study  on  the  propagation  of 
shear  waves  in  single  crystal  Cu-Al-Ni  shape  memory  alloy  has  been  done  by 
Escobar  and  Clifton  (1993).  Phase  transition  shocks  are  not  observed  directly 
due  to  their  low  propagation  speed.  Instead,  their  presence  is  inferred  from 
the  measurements  of  the  elastic  waves  at  the  rear  end  of  the  specimen.  An 
analytical  attempt  to  model  these  experiments  is  presented  in  Abeyaratne  and 
Knowles  (1997).  In  this  work  experiments  will  be  conducted  on  polycrystalline 
NiTi  SMAs. 

Classical  rate-independent  plasticity  theory  is  not  sufficient  to  describe  the 
behavior  of  SMA  materials.  While  it  is  still  capable  of  partially  predicting  the 
shape  memory  effect  (without  capturing  the  strain  recovery  upon  heating), 
it  cannot  model  the  pseudoelastic  response.  However,  for  rate  independent 
models  of  SMAs  both  theoretical  and  experimental  developments  of  dynamic 
elasto-plasticity  can  be  used  for  guidance.  Theoretical  developments  on  elasto- 
plastic  wave  propagation  in  long  rods  dates  back  to  the  works  of  Von  Kar- 
man  (1942),  Rakhmatulin  (1945)  and  Taylor  (1958).  Extensive  experiments 
on  elasto-plastic  wave  propagation  have  been  carried  out  by  Bell  (1962);  Chid- 
dister  and  Malvern  (1963);  Kolsky  (1949);  Clifton  and  Bodner  (1966),  Bodner 
and  Clifton  (1967)  using  a  split-Hopkinson  bar  apparatus.  The  split  bar  tech¬ 
nique  itself  was  introduced  by  Kolsky  (1949).  The  reader  is  referred  to  classical 
texts  on  wave  propagation  such  as  (Kolsky,  1963;  Graff,  1975)  for  additional 
information.  In  recent  years  there  have  been  extensions  to  the  Hopkinson  tech¬ 
nique  (Nemat-Nasser  et  ah,  1991)  that  allow  for  dynamic  test  recovery  in  both 
tension  and  compression.  The  basic  split-Hopkinson  technique  will  be  used  in 
this  work  to  conduct  the  dynamic  experiments  on  polycrystalline  NiTi  SMA 
rods. 


The  main  focus  of  this  paper  is  the  study  of  the  one-dimensional  dynamic 
problem  of  loading  an  SMA  rod  under  conditions  of  pseudoelasticity  and  de- 
twinning.  Both  computational  and  experimental  results  are  obtained.  Based  on 
experimental  observations  the  rate  independent  constitutiv'e  model  (Lagoudas 
et  ah,  1996)  is  selected.  The  complex  nature  of  most  constitutive  models  for 
SMA  materials  makes  direct  integration  of  even  the  simplest  uniaxial  transient 
initial  boundary  value  problems  (IB VP)  very  complicated.  Closed  form  solu¬ 
tions  can  usually  be  obtained  for  simple  boundary  conditions,  e.g.  impact  step 
loading  (Chen  and  Lagoudas,  2000)  or  by  simplifying  the  constitutive  model 
so  that  the  stress  can  be  obtained  as  an  explicit  function  of  strain  (Bekker 
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et  al  2002).  Numerical  solutions  of  the  impact  loading  problem  have  been 
obtained  by  (Jimenez- Vicktory,  1999;  Bekker  et  al.,  2002)  by  mainly  using 
the  Lax-Friedrichs  finite  difference  scheme.  This  FD  scheme  has  been  found 
to  produce  a  considerable  amount  of  numerical  dissipation  which  makes  the 
distinction  between  a  self-contained  nonlinear  shock  and  a  rarefaction  wave 
difficult.  In  this  paper  numerical  simulations  of  step  and  pulse  shock  loading 
both  for  stress  induced  phase  transformation  and  detwinning  of  martensite  are 
performed  using  the  FEM  method.  An  adaptive  meshing  technique  based  on 
the  ZZ  error  estimator  (Zienkiewicz,  1987)  is  utilized  in  order  to  improve  the 
accuracy  of  the  method  and  decrease  computational  time.  Comparisons  with 
analytical  solutions  are  made  whenever  such  solutions  are  available.  Based  on 
the  simulation  results,  the  energy  dissipation  of  SMA  rods  for  pulse  loads  are 
discussed. 

An  experiment  on  the  wave  propagation  in  a  SMA  rod  is  also  performed  in 
a  split-Hopkinson  bar  apparatus.  A  nearly  equiatomic  NiTi  SMA  specimen 
instrumented  with  strain  gauges  is  tested  under  detwinning  conditions  for  a 
pulse  impact  load.  Separate  tests  in  a  standard  uniaxial  MTS  test  frame  are 
performed  to  establish  its  quasi-static  response.  The  results  of  the  Hopkinson 
bar  experiment  are  used  to  extract  the  dynamic  stress-strain  relationship  due 
to  detwinning.  The  adaptive  FEM  technique  is  used  to  simulate  the  propaga¬ 
tion  of  stress  waves  in  the  dynamic  experiment. 

The  paper  begins  with  a  brief  overview  in  Section  2  of  the  field  equations 
and  boundary  conditions  and  constitutive  model  defining  the  problem.  The 
implementation  of  the  FEM  for  the  NiTi  SMA  is  outlined  in  Section  3.1.  The 
adaptive  strategy  is  presented  in  Section  3.2.  In  order  to  verify  the  implemen¬ 
tation  of  the  adaptive  FEM  a  boundary  value  problem  with  a  step-function 
stress  boundary  condition  is  solved  in  Section  4.1.  This  specific  boundary  con¬ 
dition  allows  for  the  construction  of  analytical  solutions  which  can  be  used 
to  verify  the  numerical  solijition  methodology.  Then,  a  square  pulse  IB  VP  is 
solved  for  conditions  of  stress  induced  martensite  (Section  4.2)  and  detwinning 
(Section  4.3).  Expected  values  for  energy  dissipation  as  the  pulse  propagates 
through  the  rod  are  presented.  Section  5  describes  the  split-Hopkinson  bar 
experrment  and  discusses  the  dynamic  characterization  of  SMA  materials.^  Fi¬ 
nally,  in  Section  5.4  the  numerical  schemes  developed  in  this  paper  are  utilized 
to  simulate  the  experimental  results. 


2  Field  equations  and  constitutive  model  for  the  impact  problem 
of  SMA  rods 


A  cylindrical  SMA  rod  of  uniform  cross-section  and  length  L  is  considered.  A 
coordinate  cover  is  associated  with  the  centroidal  axis  of  the  rod  spanning  the 
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interval  0<x  <L.  The  rod  which  is  initially  stress  free  and  at  rest  is  subjected 
to  an  impact  load  at  its  left  end  (x  =  0).  The  right  end  {x  =  L)  is  assumed  to 
remain  traction  free.  The  field  eQuations,  initial  and  boundary  conditions  are 
presented  next  followed  by  a  description  of  the  thermomechanical  constitutive 
model  for  SMAs. 


2.1  Field  equations,  initial  and  boundary  conditions 


The  rod  is  assumed  to  be  long  compared  to  its  diameter  so  it  is  under  uniaxial 
stress  state  and  the  stress  <7(2;,  t)  depends  only  on  the  axial  position  and  time. 
The  axial  component  of  the  displacement  is  denoted  by  u{x,t).  Linearized 
strain  is  further  assumed  so  the  axial  component  of  the  strain  s{x,t)  is  related 
to  the  displacement  by  e(x,  t)  =  du/dx.  Finally,  the  density  of  the  material  p 
is  assumed  constant.  The  local  form  of  the  balance  of  linear  momentum  and 
energj-  then  read  (Graff,  1975;  Malvern,  1977): 


d'^u 


da 

dx 

dx 


(1) 

(2) 


where  U  is  the  internal  energy'  per  unit  mass  and  q{x,  t)  is  the  heat  flux. 

The  timescale  of  the  impact  problem  is  on  the  order  of  micro-  to  millisec¬ 
onds.  The  physically  meaningful  IBVP  is  an  adiabatic  one  because  such  time- 
intervals  are  too  short  for  heat  conduction  to  take  place  as  well  as  for  convec¬ 
tion  to  remove  heat  through  the  surface  of  the  rod.  In  the  adiabatic  approxi¬ 
mation,  therefore,  the  heat  conduction  term  q  in  (2)  can  be  neglected  so  the 
balance  of  energy  in  conjunction  with  (1)  yields 

dU  d'^u 
^  dt  ^  dxdt 

Equation  (1)  and  (3)  involve  the  field  variables  u,  a  and  U.  Through  appro¬ 
priate  constitutive  assumptions  to  be  discussed  in  the  following  section  only 
u{x,  t)  and  the  temperature  T{x,  t)  will  become  the  independent  variables. 

For  the  field  variables  the  following  initial  and  boundary  conditions  are  as¬ 
sumed: 
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u\t=o  =  0,  ■^l«=o  =  0,  T\t=o  —  Tr 


(4) 


(^\x=o  =  a\x=L  =  0  (5) 

The  initial  conditions  indicate  that  the  rod  is  at  rest  and  its  temperature  is 
equal  to  the  ambient  temperature  Tr.  The  boundary  conditions  specify  the 
traction  ao{t)  applied  ^  to  the  left  end  of  the  rod.  The  right  end  is  kept  traction 
free. 


2.2  Thermomechanical  constitutive  model  for  poly  crystalline  SMAs 


The  field  equations  (1),  (3)  and  initial  and  boundary  conditions  (4),  (5)  alone 
are  not  sufficient  to  form  a  complete  IBVP.  A  thermomechanical  constitutive 
model  that  captures  the  key  characteristics  of  pseudoelasticity  and  detwinning 
of  the  SMA  response  is  needed. 


2.2.1  Stress  induced  martensite 

The  constitutive  model  used  is  formulated  in  terms  of  the  Gibbs  free  energy 
G  and  employs  the  volume  fraction  of  detwinned  martensite  ^  formed  from 
austenite  as  an  internal  variable  (Lagoudas  et  ah,  1996).  The  specific  form  of 
G  in  the  one  dimensional  case  is: 

G-=G{a,T,0  =  - 

+c  ((T  -  Tr)  -  Tin  (^))  -  5oT  +  uo  +  /(O 
and  it  is  linked  to  the  internal  energy  W  by  a  Legendre  transformation: 


14  =  G  +  TiS  + 

The  definition  of  G  includes  the  inelastic  transformation  strain  £*  associated 
with  the  phase  transformation.  The  function  / (^)  is  taken  to  be  a  quadratic 
polynomial  in  ^  and  is  responsible  for  the  transformation  hardening. 

^  There  is  no  continuity  requirement  on  <To(t)  i.e.  impact  loads  are  allowed 


(6) 


(7) 
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(8) 


m)  = 


\pb^^e+{^^i+^2% 

+  (Ail  -  A^2)^, 


e>o 

e<o 


where  material  constants  pb^,  pb^ ,  in  and  p2  define  the  transformation  sur¬ 
faces  and  the  hardening  during  the  forward  and  reverse  transitions  (Qidwai 
and  Lagoudas,  2000b).  In  the  above  ^  >  0  denotes  the  forward  transformation 
and  ^  <  0  the  reverse.  The  remaining  material  properties  in  (6)  are  the  effec¬ 
tive  compliance  S,  effective  thermal  expansion  coefficient  a,  effective  specific 
heat  c,  effective  specific  entropy  at  the  reference  state  Sq  and  effective  specific 
internal  energy  at  the  reference  state  uq  for  the  SMAs  which  is  composed  of 
a  mixture  of  austenite  and  martensite.  They  are  approximated  by  the  follow¬ 
ing  averaging  expressions,  which  are  good  approximations  for  polycrjstalline 
SMAs  with  random  orientation  distributions  for  grains  (Boyd  and  Lagoudas, 
1994): 


5  - 

=  +  CA5,  AS  : 

=  S^^  ■ 

-S^ 

a  — 

«(0 

=  Oi^  ^Aoi,  Aoi : 

- 

c  = 

c(0 

=  c'^  +  ^Ac,  Ac  :■ 

=  - 

(9) 

So(0 

=Z  Sq  A  ^ASq,  Asq 

._  M  . 
.—  Sq 

-  9^^ 

-uo  = 

^io(0 

=  u^  +  ^Auo,  Auo 

Uq 

- 

Uq 

Quantities  wuth  subscript  A  denote  the  appropriate  material  constant  for  the 
austenite  phase  and  those  with  subscript  M  for  the  martensite  phase.  Follow¬ 
ing  a  standard  thermodynamic  procedure  the  following  constitutive  relations 
are  obtained: 

dG 
dT 


dG 

7r  =  -p^ 


(10) 

(11) 

(12) 


where  s  is  the  entropy  and  tt  is  the  driving  force  for  the  transformation.  Using 
(11)  the  following  constitutive  relation  is  obtained: 

a  =  E{0{e-a{0{T-TR)-e^)  (13) 


where  E{^)  =  l/S{()  is  the  effective  elastic  modulus.  The  evolution  of  the 
inelastic  variable  C  is  given  by  a  consistency  condition  derived  from  a  trans¬ 
formation  criterion  (Lagoudas  et  ah,  1996).  The  evolution  of  e*  follows  that  of 


8 


(14) 


^  and  in  the  one  dimensional  case  can  be  integrated  explicitly  to  yield; 

£*  =  Hsgn{a)^ 


Here  H  is  a.  positive  material  constant  corresponding  to  the  maximum  trans¬ 
formation  strain.  The  principle  of  maximum  transformation  dissipation  in 
conjunction  with  the  second  law  of  thermodynamics  leads  to  the  following 
transformation  surface; 


TT  =  iT*  (15) 

where  y  =  -lpAso(.4”-'  -  M“)  -  ipAs„(M"  -  M”/  -  A’/  +  A-’).  The  +y 
at  the  right  hand  side  stands  for  the  forward  (A  ^  M)  transformation  surface 
and  -Y*  for  the  reverse  (M  -)•  A)  transformation  surface.  The  start  and 
finishing  temperature  for  the  forwards  transformation  are  denoted  by  A  and 
A°^  and  the  start  and  finishing  temperatures  for  the  reverse  transformation 
are  denoted  by  and  ,  respectively. 

For  detailed  description  of  the  transformation  rule  and  conditions  for  the  for¬ 
ward  and  reverse  phase  transformation  the  reader  is  referred  to  the  original 
paper  by  Lagoudas  et  al.  (1996).  The  next  section  describes  how  detwinning 
is  incorporated  into  the  constitutive  model. 

2.3  Detwinning  of  martensite 

The  detwinning  deformation  will  be  accounted  for  by  adapting  the  constitu¬ 
tive  model.  The  material  constants  for  twinned  and  detwinned  martensite  are 
the  same.  Consequently,  the  initial  response  and  the  response  after  the  com¬ 
pletion  of  detwinning  will  both  be  elastic  with  the  slope  being  the  modulus  of 
elasticity  of  martensite  .  The  deformation  is  irreversible  upon  unloading 
which,  consequently,  will  also  be  elastic. 

The  material  constants  in  the  constitutive  model  can  be  reinterpreted,  replac¬ 
ing  the  ones  for  the  austenitic  phase  with  the  ones  for  martensite.  This  will 
ensure  the  same  elastic  response  prior  the  onset  of  detwinning  and  after  its 
completion.  The  internal  variable  ^  should  be  interpreted  as  the  volume  frac¬ 
tion  of  detwinned  with  respect  to  self-accommodated  martensite  and  H  is  the 
maximum  inelastic  strain.  From  equation  (15)  the  transformation  surface  will 
have  the  following  simple  form; 


crH- 


(16) 
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The  hardening  function  in  this  case  may  be  expressed  as  follows: 

=  (17) 

where  Y‘^  =  a^H  and  pb^  =  cr/H  -  Y^.  For  convenience,  the  critical  stress 
level  Os  for  the  onset  and  Of  for  the  completion  of  the  detwinning  deformation 
are  introduced  as  material  constants.  Note  that  for  the  detwinning  case  ^  can 
only  be  positive  since  the  unloading  is  entirely  elastic.  This  adaptation  of  the 
model  allows  for  the  modelling  of  detwinning  deformations  when  no  stress 
induced  martensite  is  being  produced. 


2.4  Isentropic  approximation 


The  adiabatic  heat  equation  can  be  simplified  in  order  to  facilitate  the  numer¬ 
ical  treatment  of  the  impact  problem.  Using  the  Legendre  transformation  (7) 
the  internal  energy  can  be  eliminated  from  equation  (3). 


(18) 


Further,  upon  combining  (6)  and  (10)  an  explicit  expression  for  the  entropy 
is  obtained 


s  =  OiO ! p  -l-  cln(T/Ti?)  -h  Aso^  +  Sq  (19) 

On  substituting  (19)  into  (18)  the  balance  of  energy  becomes: 

pc^  =  {ao  -f-  pAsoC)  + 

According  to  Cory  (1985)  and  McNichols  (1987)  tt  <  pAsoT  for  most  SMAs. 
For  NiTi  the  precise  values  yield  7r/pAsoT  ^  0.013  so  equation  (20)  can  be 
approximated  by 


pc^  =  (aa  -f  pAsoO  (21) 

which  is  equivalent  to  the  isentropic  condition  St  =  0.  The  heat  capacity  c  can 
be  assumed  constant  for  the  two  phases  (i.e.  =  c^^).  Then  equation  (21) 

can  be  integrated  directly,  yielding: 
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(22) 


Consequently,  the  differential  equation  (3)  is  replaced  by  the  algebraic  equa¬ 
tion  (22).  The  impact  problem  then  reduces  to  solving  the  balance  of  linear 
momentum  (1)  for  the  only  field  variable  u{x,t).  The  remaining  field  vari¬ 
ables  a  and  T  are  coupled  with  the  strain  e  and  the  internal  variable  of  the 
constitutive  model  ^  by  equations  (13)  and  (22). 


2.5  Tangent  moduli 


A  nonlinear  displacement-based  FEM  solver  utilizing  the  Newton-Raphson 
iteration  to  resolve  the  nonlinearity  requires  partial  derivatives  of  the  stress 
with  respect  to  an  increment  of  the  strain.  An  increment  in  the  strain  causes 
increments  in  both  stress  (equation  (13))  and  temperature  (equation  (22)): 

da  _  da  da  dT 

In  order  to  find  the  total  derivative  ^  a  closed  form  expression  for  ^  is 
needed.  This  is  done  by  differentiating  equations  (13)  and  (22)  with  respect 
to  the  strain  and  combining  the  result  to  obtain: 


ot^  +  {o^oc  -f  pAso)|^j  / 

(24) 

Second  order  approximations  for  the  partial  derivatives  If  and  ^  are 

developed  in  (Qidwai  and  Lagoudas,  2000a)  and  thus  all  the  quantities  in  (23) 
can  be  computed  numerically. 


3  Numerical  implementation 


The  numerical  techniques  used  to  implement  the  constitutive  lav  s  are  de¬ 
scribed  first.  For  given  strain  increment  As  and  temperature  increment  AT 
the  stress  a  given  by  equation  (13)  is  computed  with  the  help  of  the  cutting 
plane  return-mapping  algorithm  described  in  (Qidwai  and  Lagoudas,  2000a). 
A  displacement  based  FEM  provides  strain  increments.  In  the  impact  problem 
both  stress  and  temperature  depend  on  the  strain  increment  Ae,  that  is  for 
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given  strain  both  (13)  and  (22)  have  to  be  satisfied  simultaneously.  This  is 
done  via  an  iterative  process.  The  process  starts  with  given  values  cr(°^ 
T^^'>  for  strain,  stress  and  temperature  which  satisfy  (13)  and  (22).  Given  a 
strain  increment  Ae  the  pair  (a,T)  corresponding  to  strain  e  =  +  Ae  is 

found  through  the  iteration: 

^(n+i)  =,E(e-a  -  Tr)  -  (25) 

y(n+l)  ^2^^g-^(oa("+‘)+pAsoCC"+‘))  (26) 


The  first  equation  (25)  uses  the  return-mapping  algorithm  to  compute  a  new 
value  for  the  stress  based  on  the  old  temperature  The  second 

equation  (26)  attempts  to  enforce  the  isentropic  heat  equation  by  comput¬ 
ing  a  corrected  temperature  The  process  is  terminated  when  there  is 

no  further  progress,  i.e.  when  I  and  both  become 

smaller  than  certain  tolerance.  The  algorithm  showed  linear  convergence  in 
the  test  cases,  however  a  detailed  theoretical  study  is  required  to  establish  its 
properties. 


3.1  FEM  procedure 


A  standard  semi-discrete  Galerkin  approximation  is  used  to  generate  the  weak 
form  of  the  problem.  In  this  paper  only  linear  elements  will  be  used.  Let 
P^([0,L])  C  H^{[0,L])  be  the  set  of  piecewise  linear  functions  over  each  ele¬ 
ment  and  {V'ilili  be  the  usual  basis  of  P^([0,  L]).  The  weak  form  of  (1)  then 
reads: 


Find  u^{x,t)  =  Eili  Ui{t)'ipi{x)  such  that  for  Vu'"  €  P^{[0,  L]): 


rL  52^/1  r 

>’fo 


L  Qy^ 
a—ir^dx 


dx 


—av 


h 


x=0 


(27) 


As  usual  the  number  of  nodes  is  N  (i.e.  N—\  elements)  and  the  nodal  values  for 
the  displacement  are  denoted  by  Ui{t).  Whenever  appropriate,  vector  notation 
will  be  used,  that  is  U  =  (C/j, ..., Piv)^  Problem  (27)  is  reduced  to  a  second 
order  nonlinear  system  of  ODEs: 

MU  =  F(U)  (28) 


where  M  is  the  mass  matrix  and  F^(t)(U)  is  the  forcing  term.  The  subscript 
^{t)  stands  to  indicate  that  F^(t)(U)  does  not  depend  on  the  displacement 
only  but  on  the  whole  loading  history.  However,  for  any  given  loading  history 
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the  stress  and  hence  F|(t)(U)  can  be  viewed  as  well  defined  single  valued 
functions.  Thus,  without  loss  of  generality  the  subscript  ^{t)  will  be  dropped 
in  the  discussion  that  follows.  The  mass  matrix  and  load  vector  are  given  by: 


(29) 

(30) 


It  is  also  useful  to  introduce  the  forcing  term  F(U)  due  to  inelastic  strains 
and  the  stifllness  matrix  K(U)  which  are  given  by^  : 


F,{V)  =  E{(.)  [e'(0  +  a«)(r  -  li)]  ^dx 


(31) 

(32) 


Note  that  the  decomposition  F(U)  =  F(U)  -  K(U)U  holds  and  (28)  can  be 
rewritten  as: 

MU  +  K(U)U  =  F(U)  (33) 


The  time  integration  in  (28)  (or  (33))  is  done  by  the  backward  difference 
method,  a  member  of  the  Newmark  family  (Newmark,  1959;  Reddy,  1993). 
For  t  —  ts  the  Newmark  scheme  is  defined  by^  : 


Urt,  =  U,  +  tU.  +  ir^((l  -  7)U,  +  7Urt,)  (34) 

Uji+1  =  Us  +  r((l  -  q;)Us  +  aUs+i)  (35) 

The  backward  difference  method  is  obtained  by  setting  a  =  |  and  7  =  2.  It  is 
easy  to  show  (see  e.g.  (Reddy,  1993))  that  the  above  difference  equations  lead 
to  the  following  system  of  nonlinear  algebraic  equations  for  [/s+i: 

=  F(U.+,)  +  G.  (36) 

or,  equivalently,  to 


.4M  +  K(U,+,))u,+,  =  F(Urt,)  +  G,  (37) 

7t2  j 


2  Similarly,  a  more  precise  notation  for  K  and  F  would  be  K^(q(U)  and  F^(t)(U), 
respectively. 

3  The  usual  notation  Us  :=  U(ts)  is  used 
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where  Gj  =  M  (:^Us  +  :^Us  +  ^Ui)-  The  nonlinear  problem  (36)  is  solved 
by  linearizing  the  right-hand  side 


N 


Fi{V  +  AV)c,Fi{V)  +  Y: 

i=i 


dFijU) 

dUj 


AUj 


and  using  the  chain  rule  to  obtain: 


Lij{U) 


dFijV) 

dUj 


da  dipi  _  [^  da  d'lpj  dipi 


(38) 


The  solution  Us+i  is  found  through  a  Newton-Raphson  iterative  process.  Set 
the  initial  guess  to  and  forn  =  1, 2 . . .  until  convergence  compute: 


-  L(U<’;>,)  j  (F(uS”>i)  -  L(uW,)u‘?i  +  G.)  (39) 

The  cutting  plane  method  (Qidwai  and  Lagoudas,  2000a)  which  is  used  to  re¬ 
solve  the  nonlinear  behavior  of  the  material  also  provides  second  order  numer¬ 
ical  approximation  for  the  derivative  dajde  which  results  in  a  quasi-Newton 
algorithm.  Since  the  Newton  algorithm  is  only  locally  convergent  in  the  cases 
when  it  diverges  the  simple  iteration  was  applied  to  (37).  Again,  set  =  Uj 
and  for  n  =  1, 2 . . .  until  convergence  compute: 


+  K{U<”>0U<”',  j  (F(uS,)  +  G.)  (40) 

In  all  numerical  examples  tested  the  later  iteration  demonstrated  global  linear 
convergence. 


3.2  Adaptive  mesh  refinement 


Let  a^  be  the  stress  at  the  completion  of  the  Newton  iterations  for  given  time 
step  n,  i.e.  t  —  tn-  For  linear  elements  a*^  is  a  piecewise  constant  function.  Let 
be  the  continuous,  piecewise  linear  function  in  [0,  L]  which  assumes  the 
averaged  value  of  a^  at  each  nodal  point.  The  error  indicator  r]^{e)  is  defined 
locally  over  each  element  e  by  (Zienkiewicz,  1987): 


ria{&)  = 


(41) 


where  H-Hog  is  the  L2  norm.  An  element  e  is  refined  if 

r]a{e)famax  >  TOLl 


(42) 
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where  cXmax  is  the  absolute  value  of  the  maximum  attainable  stress  in  the  rod, 
which  for  impact  problems  is  known  in  advance.  Two  neighboring  elements  Cj 
and  Ci+i  are  merged  into  one  if 

r}Aei)l^max  <  TOL2,  Vaiei+i)lamax  <  TOL2  (43) 

Two  aspects  of  the  actual  implementation  details  of  the  FE  analysis  should  be 
emphasized.  The  linear  system  (36)  (or  (37))  is  tridiagonal  and  poses  no  com¬ 
putational  problems.  Secondly,  the  most  time-consuming  parts  of  the  FE  pro¬ 
cedure  are  the  assembly  of  the  stiffness  matrix  at  each  Newton  step  (because 
of  the  nonlinear  dependance  of  the  stiffness  on  the  strain)  and  the  assembly  of 
the  force  vector.  They  require  the  execution  of  the  stress  update  procedure  via 
the  return-mapping  algorithm  which  is  a  computationally  expensive  operation 
and  is  preformed  once  for  each  element  at  each  Newton  step. 

Clearly  a  global  uniform  h-refinement  strategy’  used  to  achieve  satisfactory 
spatial  discretization  will  impose  severe  restrictions  on  the  problem  size  due 
to  the  assembly  time  issues.  In  order  to  avoid  this  the  local  criterion  (42)  is 
applied  to  each  element  at  the  completion  of  the  Newton  iteration  to  refine  or 
coarsen  the  mesh.  If  there  is  no  further  need  to  refine  the  mesh  the  algorithm 
proceeds  to  the  next  time  step.  It  was  found  that  this  approach  works  very 
well  for  the  class  of  SMA  hysteretic  materials  under  consideration. 


4  Numerical  Examples 


The  implementation  of  the  FEM  was  tested  in  three  different  numerical  exam¬ 
ples.  The  step  loading  problem  under  conditions  of  pseudoelasticity  (T  >  A°^) 
presented  in  the  next  section  is  used  to  compare  the  numerical  solution  to  ex¬ 
isting  analytical  solutions  (Chen  and  Lagoudas,  2000;  Bekker  et  ah,  2002). 
It  is  also  used  to  demonstrate  the  capabilities  of  the  adaptive  mesh  refine¬ 
ment  strategy.  Secondly,  a  problem  with  pulse  boundary  conditions  is  solved, 
again  under  pseudoelastic  conditions.  The  third  problem  also  features  a  pulse 
boundary  condition  but  at  a  lower  temperature  (T  <  M°®)  so  only  detwinning 
of  martensite  is  involved. 

The  material  properties  (Table  1)  for  all  model  problems  are  taken  from  (Qid- 
wai  and  Lagoudas,  2000a)  and  represent  generic  NiTi  SMA  properties.  In  ad¬ 
dition  to  that  for  all  numerical  simulations  the  length  of  the  rod  was  taken  to 
be  0.5m.  All  calculations  were  performed  on  a  933  Mhz  PHI  machine  running 
Windows  NT. 
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Table  1 


Material  parameters  used  in  the  SK 

:A  model 

Material  constant 

Value 

Material  constant 

Value 

70  X  10®  Pa 

dcr 

dT 

7.0  X  10®  Po/(m®A) 

30  X  10®  Pa 

275  °K 

22  X  IQ-^/K 

291  °K 

10  X  10-^ /K 

A"® 

295  °K 

H 

0.05 

A°1 

315  °K 

4.1  Step  loading  problem 


The  fixed  impact  stress  initial-boundary  value  problem  ^  is  defined  by  setting 
the  boundary  condition  to  be  the  step  function: 


f  0  for  i  <  0 
I  ao  for  t  >  0 


(44) 


The  strain  level  eo  which  causes  the  constant  impact  stress  gq  can  be  found 
from  equation  (13).  This  particular  boundary  condition  is  chosen  because  it  is 
a  natural  starting  point  for  nonlinear  hyperbolic  equations  and  because  there 
are  existing  analytical  solutions  for  it. 


4.1.1  Analytical  solutions  to  the  step  loading  problem 

The  structure  of  the  solution  depends  strongly  on  the  impact  stress  tJo.  Let  the 
pair  (se/,  cTei)  be  the  point  on  the  hysteresis  curve  that  corresponds  to  the  start 
of  the  phase  transformation.  In  this  example  cq  it  is  taken  to  be  sufficiently 
high  so  that  full  phase  transformation  transformation  has  occurred.  It  is  also 
required  that  the  value  of  ao  be  high  enough,  so  that  the  graph  of  of  the  stress 
strain  relationship  of  the  SMA  is  below  the  line  connecting  the  points  {cei,  <Jei) 
and  (£0)<7o)  (see  Figure  1). 

Following  Chen  and  Lagoudas  (2000);  Bekker  et  al.  (2002)  it  can  be  shown 
that  for  material  with  initial  linear  stress-strain  relationship  prior  to  the  onset 
of  phase  transformation  the  solution  has  the  following  two-shock  structure: 


^  When  the  same  initial  boundary  value  problem  is  reformulated  as  an  initial  prob¬ 
lem  on  an  infinite  domain  with  the  initial  condition  being  a  step  function  it  is  usually 
referred  to  as  the  Riemann  problem. 
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Fig.  1.  Schematic  of  the  loading  portion  of  a  stress-strain  relationship  and  the  critical 
points  defining  the  solution  to  the  problem. 


a{x,t)  = 


(To  for  0  <  x/t  <  Vph 
<  Cel  for  Vph  <xft<  Vel 
0  for  Vel  <  x/t 


T{x,t) 


/ 

Tq  for  0  <  x/t  <  Vph 
<  Tel  for  Vph  <  x/t  <  Vel 
0  for  Vel  <  x/t 


(45) 


(46) 


where  To  is  the  temperature  corresponding  to  the  impact  stress  cto  and  Tei 
is  the  temperature  just  prior  to  the  onset  of  the  phase  transformation.  The 
faster  shock  is  a  linear  thermoelastic  elastic  shock  and  has  velocity 


K,  =  ./^  (47) 

V  P^el 

This  shock  is  due  to  the  shock  type  of  the  boundary  condition  and  the  initial 
linear  stress-strain  response.  The  second,  slower  shock,  is  a  transformation 
shock  which  travels  with  velocity 


I  (Tq  —  agl 

p{eo  —  Set) 


(48) 


This  shock  occurs  not  only  because  of  the  boundary  condition  but  also  be¬ 
cause  of  the  convex-dowm  nature  of  the  stress-strain  relationship  for  e  >  £ei. 
Higher  stress  levels  travel  with  higher  velocity  than  lower  stress  levels  which 


17 


make  the  shock  self  sustained  and  independent  of  the  boundary  condition  (see 
(Godlewsky  and  Raviart,  1996,  pg.  83-97)  for  a  general  discussion  as  well  as 
(Chen  and  Lagoudas,  2000;  Bekker  et  al,  2002)  for  solutions  specific  to  SMA 
materials).  The  phase  transformation  shock  specifies  the  point  of  abrupt  phase 
transition.  For  material  points  with  x  <  Vpht  the  material  is  in  the  martensitic 
phase  and  the  region  x  >  Vpht  is  still  in  the  austenitic  phase. 

Note  that  the  adiabatic  heat  equation  (22)  does  not  provide  for  a  completely 
linear  initial  response.  However,  prior  to  the  onset  of  phase  transformation, 
^  =  0  and  the  heat  equation  (22)  can  be  linearized  as  follows: 

T  =  Tn{l-^a)  +  o({^f)  (49) 

pc  \  pc  J 

By  neglecting  the  higher  order  terms  in  (49)  the  remaining  linear  part  can 
be  substituted  in  (13)  to  obtain  a  completely  linear  adiabatic  stress-strain 
response.  The  linear  approximation  in  (49)  is  justified  in  the  thermoelastic 
range  before  commencement  of  phase  transformation  because  ^  w  10“^.  If 
equation  (22)  is  not  linearized  the  elastic  shock  will  be  replaced  by  a  continuous 
function  with  very  high  gradient.  The  velocity  of  the  points  on  the  graph  of 
this  function  will  deviate  from  the  velocity  of  the  elastic  shock  by  «  10“°. 


4.1.2  Numerical  results  for  the  step  loading  problem 

For  all  numerical  simulations  the  impact  stress  level  is  ao  =  — 400MPa  cor¬ 
responding  to  impact  strain  of  e  =  0.0635.  The  reference  temperature  is 
Tr  =  320  °K.  The  FEM  solver  was  set  to  use  the  backward  difference  time 
integration  scheme  and  the  Newton-Raphson  method  to  solve  the  nonlinear 
system  (36).  The  Newton-Raphson  iteration  showed  quadratic  convergence  at 
all  time  steps  except  for  the  first  few  ones  when  the  shock  were  forming.  In 
the  cases  when  it  was  diverging  the  alternative  direct  iteration  (40)  approach 
was  used. 

Significant  computational  savings  can  be  obtained  if  isothermal  instead  of  adi¬ 
abatic  conditions  are  assumed.  In  an  isothermal  problem  the  temperature  is 
held  constant  T  =  Tr  and  the  balance  of  energy  (2)  is  not  considered.  Thus 
the  quasi-static  hysteresis  of  the  material  is  used  instead  of  solving  equations 
(13)  and  (22).  For  a  NiTi  SMA  with  the  material  data  from  Table  1  the  dif¬ 
ference  between  the  adiabatic  and  isothermal  hysteresis  is  shown  in  Figure  2. 
The  shape  of  the  hysteresis  is  the  same  and  the  differences  in  the  transforma¬ 
tion  portion  will  not  affect  the  structure  of  the  solution  provided  that  ao  is 
well  above  the  stress  level  required  to  finish  the  transformation.  Consequently, 
no  matter  whether  isothermal  or  adiabatic  conditions  are  assumed  the  shock 
speeds  Vph  and  Ki  will  only  depend  on  the  values  for  Sei,  Oei,  Eq,  uq.  From 
a  computational  point  of  view  this  simplification  avoids  the  iteration  process 
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Strain 

Fig.  2.  An  adiabatic  and  isothermal  path  for  the  material  data  in  Table  1  at 
T  =  320  °K.  Under  adiabatic  conditions  higher  stress  levels  are  required  to  complete 
the  phase  transformation  compared  to  isothermal  hysteresis  loops. 

(25), (26)  (typically  6-7  iterations)  which  results  in  a  significant  reduction  in 
computational  time.  While  the  structure  of  the  solution  is  not  compromised 
very  fine  spatial  meshes  can  be  explored  for  the  purposes  of  comparing  ana¬ 
lytical  and  numerical  solutions. 

For  the  isothermal  hysteresis  (Figure  2)  an  impact  stress  of  ao  =  400 Af Pa  is 
sufficient  for  the  full  completion  of  the  phase  transformation  under  isothermal 
conditions.  The  onset  of  phase  transformation  begins  at  a^i  =  —19bMPa  for 
a  strain  Sei  =  2.78  x  10~^.  Given  this,  the  speed  of  the  two  shocks  (48)  and 
(47)  are  found  to  be: 


Vj,h  =  723m/s  (50) 

Kj  =  3294m/s  (51) 

Based  on  the  first  few  numerical  results  (Figure  3)  and  (Figure  4)  several  ob¬ 
servations  can  be  made.  First,  all  numerical  solutions  have  the  expected  two 
shocks  -  one  elastic  and  another  corresponding  to  the  phase  transformation. 
Fixed  meshes  with  coarse  spatial  discretizations  have  oscillations  close  to  the 
phase  shock  location.  A  comparison  of  the  two  meshes  in  Figure  3,  both  for 
a  fixed  time-step  of  r  =  O.l^s  at  time  t  =  30/is  shows  that  oscillations  can 
be  eliminated  by  refining  the  mesh.  Secondly,  the  backward  difference  scheme 
which  was  used  in  these  computations,  introduced  numerical  dissipation  which 
is  most  pronounced  at  the  elastic  shock.  Several  other  members  of  the  New- 
mark  family  were  tested.  Explicit  methods  as  well  as  the  constant  acceleration 
scheme  were  found  to  be  unconditionally  unstable  producing  highly  oscillatory 
solutions  that  were  diverging  with  time.  Of  those  methods  that  were  able  to 
converge  the  backward  difference  was  found  to  dampen  the  high  frequency 
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Fig.  3.  Stress  profile  at  30/is  for  a  fixed  mesh  with  500  (a)  and  2000  elements 
(b).  Numerical  oscillations  are  eliminated  for  the  finer  spatial  discretization.  The 
position  of  the  elastic  shock  is  marked  by  a  dashed  line. 


(a)  (b) 

Fig.  4.  Stress  profile  at  30/is  for  an  adaptive  mesh  with  two  different  time  steps. 
The  linear  shock  is  smeared  for  a  coarse  time  step  r  =  0.1/^s  (a).  It  is  much  sharper 
when  a  finer  step  of  r  =  0.001/us  (b)  is  used.  Mesh  nodes  are  marked  with  black 
squares  and  the  thin  line  at  the  top  shows  the  density  of  elements. 

oscillations  (Figure  3(a))  in  the  most  efficient  manner  and  was  subsequently 
chosen  for  all  future  computations.  The  numerical  dissipation  can  be  decreased 
by  appropriately  decreasing  the  time  step.  The  quasi-Newton  method  used  to 
solve  the  nonlinear  system  (36)  showed  quadratic  convergence  at  all  time  steps 
but  the  first  few  ones  when  the  shock  were  forming.  In  that  case  the  alternative 
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direct  iteration  (40)  approach  was  used. 

Quantitatively  the  results  obtained  by  both  the  fixed  and  adaptive  FEM  are 
in  agreement  with  the  analytical  solution.  In  regions  away  from  the  shocks 
the  relative  difference  in  the  values  of  the  stress  for  the  numerical  and  the 
exact  solution  is  less  then  10““^.  The  accuracy  of  the  solutions  therefore  is 
determined  based  on  the  quality  of  the  numerical  solution  close  to  the  shock 
locations.  The  interval  covering  a  shock  (phase  or  elastic)  w'here  the  numerical 
values  for  the  stress  differ  from  the  exact  ones  by  more  than  1%  is  assumed 
to  be  the  range  of  uncertainty  for  the  numerical  value  of  the  shock  location. 
Consequently,  the  left  and  right  end  of  this  interval  are  assumed  to  be  bounds 
for  the  position  of  the  shock  of  the  numerical  solution. 

Based  on  this  measure  of  error,  for  a  time  step  of  r  =  0.1/rs  the  phase  shock 
is  found  to  travel  with  velocity  in  the  range  693  —  900m/ s.  The  velocity  of 
the  elastic  front  is  calculated  to  be  in  the  range  3316  ±  4207n/s.  These  results 
are  the  same  for  a  fixed  (Figure  3(b))  and  adaptive  mesh  (Figure  4(a)).  This 
indicates  that  the  adaptive  and  fixed  FEM  converge  to  the  same  solution. 

The  smearing  of  the  stress  profile  in  the  region  of  the  elastic  shock  is  due  to  the 
time-integration  scheme.  When  the  time  step  is  decreased  the  slope  becomes 
steeper  and  eventually  converges  to  the  shock.  For  an  adaptive  solution  with 
a  time  step  r  =  O.OOl/us  (the  same  computation  for  a  fixed  mesh  was  time 
prohibitive)  the  calculated  values  for  the  phase  shock  are  now  in  the  range 
723  —  733m/s  and  the  elastic  shock  is  within  the  bounds  3256  —  3366m/s 
(Figure  4(b)).  This  indicates  that  the  lower  bound  for  the  transformation 
shock  is  very  close  to  the  actual  value  (50)  and  that  the  elastic  shock  (51)  is 
virtually  in  the  middle  of  the  suggested  numerical  range.  The  relative  error 
in  the  predicted  value  for  the  phase  shock  velocity  decreases  from  24%  for 
r  =  0.1/iS  down  to  1.3%  for  r  =  0.001/iS.  The  error  in  the  elastic  shock 
speed  decreases  from  12%  to  1.1%  which  is  a  clear  indication  that  the  FEM 
algorithm  is  converging  to  the  exact  solution. 

An  inspection  of  Figure  3  reveals  that  there  are  large  regions  in  the  bar  with 
no  variation  in  the  stress.  This  is  fully  utilized  by  the  adaptive  approach. 
Figure  4(a)  shows  an  adaptive  FE  solution  with  the  same  time  step  as  the 
solution  on  Figure  3(b)  and  a  adaptive  tolerance  (see  (42))  set  to  lO"^.  This 
accuracy  is  comparable  to  the  one  of  a  fixed  mesh  with  2000  elements.  The 
maximum  number  of  elements  that  the  adaptive  mesh  contained  was  305.  The 
order  of  magnitude  fewer  number  of  elements  in  the  adaptive  meshes  induced 
a  corresponding  order  of  magnitude  decrease  in  the  computational  time. 

A  comparison  in  the  performance  of  the  fixed  and  adaptive  FE  methods  is 
given  in  Table  2.  The  time  step  is  r  =  0.01/us.  The  number  of  elements  for  the 
fixed  FEM  is  16000.  The  adaptive  solution  was  chosen  so  that  it  had  compa- 
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Table  2 


Execution  times  for  fixed  and  adaptive  meshes 


Time 

Fixed  Mesh 

Adaptive  Mesh 

Elements 

Time  (min) 

Elements 

Time  (min) 

10  ps 

16000 

56 

161 

1:12 

20  ps 

16000 

113 

199 

2:37 

- j 

40  ps 

16000 

226 

256 

6:10 

80  ps 

16000 

451 

301 

15 

rable  accuracy  with  the  one  for  the  fixed  mesh  solution.  A  comparison  of  the 
execution  times  for  the  fixed  and  adaptive  methods  shows  that  the  adaptive 
procedure  delivers  an  order  of  magnitude  improvement  in  performance. 


J1..2  Square  pulse  loading  problem  in  pseudoelastic  conditions 


A  more  realistic  initial-boundary  value  problem  is  one  for  which,  instead  of 
step  loading,  the  boundary  condition  is  a  square  pulse,  that  is 


ao{t)  =  < 


0  for  t  <  0 

(To  for  0  <t  <  tpulse 

0  for  t  ^  Ipuisc 


(52) 


where  tpuise  is  the  duration  of  the  pulse.  Due  to  the  complicated  constitutive 
response  and  boundary  conditions  there  is  no  analytical  solution  to  be  com¬ 
pared  with.  Moreover,  there  are  unresolved  questions  regarding  the  uniqueness 
of  the  weak  solution  for  times  t  >  tpuise  when  unloading  takes  place. 

The  stress  level  used  for  the  numerical  simulation  is  (Tq  =  800MPa  and  the 
initial  temperature  is  Tr  =  320  °K  >  A°^ .  The  simulation  is  done  for  adiabatic 
conditions,  utilizing  both  equations  (13)  and  (22)  to  calculate  the  adiabatic 
response  of  the  SMA.  The  stress  level  is  chosen  so  that  the  full  adiabatic 
hysteresis  loop  can  be  realized  (see  Figure  2).  The  pulse  length  is  tpuise  =  lOps 
and  the  time  step  is  t  =  0.001/rs. 

The  evolution  of  the  stress  and  temperature  in  the  rod  up  to  90/iS  is  shown 
in  Figures  5  and  7.  As  predicted  by  (45)  the  two-shock  solution  for  the  stress 
is  clearly  visible  at  the  end  of  the  pulse  load  at  t  =  10 ps  (Figure  5).  The  tem¬ 
perature  profile  (Figure  7)  also  has  two  shocks  (equation  (46)).  The  maximum 
temperature  To  =  378.8  °K  is  achieved  in  the  region  of  full  phase  transforma- 
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Fig.  5.  Stress  profile  at  different  instances  of  time  for  a  square  pulse  in  adiabatic 
loading 


0.00  0.01  0.02  0.03  0.04  0.05 

Axial  position  (meters) 

Fig.  6.  Magnified  view  near  the  left  end.  The  unloading  (lO/ts)  produces  two 
right-travelling  shock  waves  (20^s).  The  faster  unloading  wave  reflects  off  the  trans¬ 
formation  shock  (w  2lfis)  and  forms  a  left-travelling  wave  (24;is).  What  follows  is 
a  series  of  complicated  reflections  that  gradually  kill  the  initial  non-linear  shock. 

tion.  The  jump  in  the  elastic  shock  is  Tgi  -Tr  —  0.66  °K  and  for  this  reason 
it  is  not  clearly  visible  in  the  figure. 

The  most  noticeable  feature  observed  in  Figure  6  is  the  structure  of  the  un¬ 
loading  pulse.  Again  a  two  wave  shock  structure  is  seen  that  corresponds  to  the 
initial  elastic  unloading  and  the  following  reverse  transformation  M‘  ^  A  as 
can  be  seen  from  the  stress  profile  at  10  and  20/is.  Both  unloading  shocks  travel 
faster  than  the  forward  phase  transformation  shock.  When  the  faster  unload¬ 
ing  front  catches  up  with  the  forward  phase  transformation  shock  {t  «  22/2s) 
a  left-travelling  reflection  is  generated.  The  left-travelling  wave,  as  seen  for 
t  =  24//S,  partially  reflects  from  the  slower  unloading  shock  and  partially  con¬ 
tinues  [t  =  2^jis)  until  it  reflects  off  the  left  end  of  the  rod.  A  complicated 
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Fig.  7.  Temperature  profile  at  various  times.  The  jump  at  the  forward  transformation 
shock  is  To  —  —  58.2  °K.  The  elastic  shock  is  not  visible  clearly  because  of  its 

small  magnitude  of  T^i  —  Tr  =  0.66  °K. 


series  of  reflection  waves  follows.  The  first  reflection  results  in  approximately 
34%  decrease  of  the  peak  stress  level  {t  =  24/ts).  The  picture  becomes  even 
more  complicated  when  the  slower  unloading  shock  eventually  catches  up  with 
the  forward  travelling  phase  transformation  shock.  Eventually  the  peak  stress 
levels  are  reduced  to  values  below  Oei,  the  critical  stress  corresponding  to  the 
onset  of  phase  transformation.  The  temperature  profile  at  t  =  90/xs  is  hardly 
visible  because  the  material  is  entirely  in  the  elastic  range  and  the  temperature 
in  the  rod  is  very  close  to  the  reference  temperature.  The  large  amounts  of 
latent  heat  generated  during  the  initial  loading  phase  are  gradually  consumed 
in  the  reverse  transformation  as  the  stress  is  reduced  within  the  elastic  limits. 


For  pulse  loading  it  is  physically  meaningful  to  compute  the  energy  dissipation 
due  to  the  phase  transformation.  If  P{t)  is  the  work  done  by  the  external  forces 
at  the  left  end  of  the  rod  from  f  =  0  up  to  f  =  r,  /C(r)  is  the  kinetic  energy 
of  the  rod  at  time  t  =  T  and  W(r)  is  the  stored  elastic  energy  of  the  rod  then 
the  energy  dissipation  is  defined  by 

P(r)  -  (X(r)  +>V(t)) 

P(t) 


The  quantities  P,  K,  and  W  given  by 

fo''cr(0,t)v(0,t)dt 

W(r)-^fo^a(x,T)£^(x,r)dx  (54) 
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Fig.  8.  Energy  dissipation  for  a  10  /rs  square  pulse  in  adiabatic  conditions, 
can  be  easily  computed  numerically  at  each  time  step. 

The  calculations  show  (Figure  8)  that  the  dissipation  level  goes  from  40% 
at  the  end  of  the  pulse  (T  =  lO^s)  to  64%  at  T  22/is  when  the  faster 
unloading  wave  reflects  off  the  forward  travelling  transformation  wave.  The 
high  stress  levels  are  then  gradually  reduced  within  the  elastic  limits.  The 
energy  dissipation  reaches  approximately  84%  at  100  /zs,  shortly  before  the 
elastic  front  reaches  the  right  end. 


4 .3  Detwinning  Induced  by  a  pulse  load 


In  this  numerical  simulation  the  same  boundary  condition  (52)  as  in  the  pre¬ 
vious  section  is  used.  The  initial  temperature  is  set  to  Tr  =  295  °K  which 
is  in  the  detwinning  range  and  the  material  is  initially  in  the  state.  The 
stress  pulse  has  magnitude  oq  =  AOOMPa  which  is  sufficient  to  complete  the 
detwinning  and  then  obtain  the  elastic  response  of  the  martensite  phase. 

There  is  no  latent  heat  generation  during  the  detwinning  deformation.  If  it  is 
assumed  that  all  the  work  dissipated  through  inelastic  deformations  is  trans¬ 
formed  into  heat,  then  the  change  in  temperature  would  be  2  °K.  Therefore 
it  is  both  physically  and  computationally  justified  to  perform  the  simulation  in 
an  isothermal  setting.  The  loading  part  of  the  hysteresis  is  of  the  same  type  as 
the  loading  part  {A  — >  M^)  of  the  stress-strain  relationship  for  stress  induced 
martensite.  Therefore  for  the  duration  of  the  pulse  a  two-shock  structure  for 
the  stress  distribution  can  be  expected  (see  equations  (45),  (50)  and  (51)). 
This  is  observed  clearly  for  the  stress  profile  at  t  =  10/zs  in  Figure  10. 

The  unloading  is  completely  elastic  and  a  single  linear  shock  forms,  travelling 
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Fig.  9.  Stress  profiles  at  various  times  for  a  square  pulse.  The  material  is  the  de- 
twinning  range.  The  attenuation  of  the  stress  to  values  within  the  elastic  material 
response  is  clearly  visible. 
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Fig.  10.  Magnified  view  of  stress  profiles  in  the  region  close  to  the  left  end  of  the 
rod.  Direction  of  the  shock  velocities  are  indicated  by  arrows. 


at  the  speed  of  the  forward  elastic  shock  (both  the  initial  loading  response  and 
unloading  are  linear  with  the  elastic  modulus  of  martensite).  The  unloading 
shock  is  therefore  fast  enough  to  catch  up  with  the  nonlinear  shock  caused  by 
the  detwinning.  This  is  followed  by  a  series  of  reflections  between  the  left  end 
(which  is  traction  free  after  the  pulse  is  over)  and  the  forward  propagating 
detwinning  shock.  The  stress  profile  at  several  diflferent  instances  of  time  is 
presented  in  Figure  9. 


The  energy  dissipation  (Figure  11)  in  the  rod  follows  a  similar  path  as  in  the 
previous  numerical  simulation.  The  first  significant  rise  in  the  dissipation  levels 
occurs  immediately  after  unloading,  at  i  =  lO/xs.  After  the  unloading  wave 
reaches  the  forward  propagating  detwinning  front  at  i  w  18/zs  a  new  rise  in 
the  dissipation  occurs  leading  to  final  levels  of  approximately  86%.  It  should 
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Fig.  11.  Energy  dissipation  for  a  10  fis  square  pulse. 

be  noted  that  this  case  is  not  equivalent  to  the  pulse  load  in  pseudoelastic 
conditions  because  of  different  initial  stress  levels.  Another  difference  with  the 
pseudoelastic  case  is  that  the  material  is  permanently  deformed  and  in  order 
to  recover  its  shape  the  rod  has  to  be  reheated. 


5  Dynamic  Loading  Experiment  of  an  SMA  rod 


The  dynamic  response  of  a  nearly  equiatomic  NiTi  alloy  rod  is  character¬ 
ized  with  one  dimensional  wave  propagation  experiments  in  a  Hopkinson  bar 
arrangement.  The  main  feature  of  the  Hopkinson  implementation  of  the  dy¬ 
namic  experiment  is  in  the  length  of  the  specimen,  L^p  which  is  quite  long. 
This  means  that  a  steady-state  condition  is  not  reached  during  the  time  of 
the  experiment  and  one  has  to  deal  with  the  propagation  of  the  wave  in  the 
specimen  material. 


5.1  Description  of  the  Apparatus 


Hopkinson  bar  apparatus  has  become  standard  in  the  characterization  of  the 
dynamic  response  of  materials.  Detailed  descriptions  are  provided  in  many 
handbooks  and  textbooks  (Kolsky,  1963;  Graff,  1975),  and  hence  only  a  brief 
description  is  provided  here.  A  photograph  of  the  experimental  setup  is  shown 
in  Figure  12  and  a  schematic  of  the  impact  device  is  given  in  Figure  13. 

The  apparatus  consists  of  a  striker  bar,  an  input  bar  and  an  output  bar,  all  of 
diameter  d  =  15.5mm  and  all  made  of  a  4340  steel,  quenched  and  tempered 
to  a  martensitic  state.  The  yield  strength  of  these  bars  is  about  1.8  GPa  and 
they  remain  elastic  during  the  impact  experiments.  The  density  of  the  bars  is 
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Fig.  12.  Photograph  of  the  Hopkinson  bar  experimental  setup.  The  specimen  is 
visible  at  the  top-right  part  of  the  photograph 


Incident  pulse 


Fig.  13.  Geometry  and  arrangement  of  strain  gauges  in  Hopkinson  apparatus.  (Fig¬ 
ure  not  drawn  to  scale) 

p  —  TSOO/c^/m^,  the  measured  bar  wave  speed  Ct,  =  sjE^/p  =  5300m/s  and 
Eb  is  the  modulus  of  elasticity  of  the  steel  bar.  The  striker  bar  (13)  of  length  L 
is  propelled  from  an  air  gun  at  speeds  in  the  range  of  10  to  40  m/ s.  This  striker 
impacts  the  input  bar  which  is  1.7m  long.  A  one  dimensional  compression  wave 
propagates  into  both  bars.  Since  the  striker  bar  is  short,  the  reflected  tension 
pulse  arrives  at  the  striker-input  bar  interface  at  a  time  tpuUe  =  ‘^ElCt-  At 
this  point,  the  striker  comes  to  a  stop  and  is  disengaged  from  the  input  bar. 
Hence,  a  compression  pulse  of  duration  tpuise  is  propagated  down  the  length 


of  the  input  bar.  This  wave  is  coupled  into  the  specimen  which  is  in  contact 
with  the  far  end  of  the  input  bar.  Due  to  the  impedance  mismatch  between 
the  specimen  and  the  input  bar,  part  of  the  pulse  is  reflected  back  into  the 
input  bar  and  part  of  the  pulse  propagates  into  the  specimen.  A  strain  gauge 
mounted  at  about  the  middle  of  the  input  bar  is  used  to  monitor  the  incident 
compressive  pulse  and  the  reflected  tensile  pulse  propagating  in  the  input  bar. 
The  wave  propagating  through  the  specimen,  gets  coupled  into  the  output  bar, 
again  with  a  reflected  component  due  to  the  impedance  mismatch.  The  output 
bar  is  free  at  the  far  end  and  so  a  tensile  pulse  reflects  from  the  far  end  of  the 
output  bar  and  is  unable  to  transmit  into  the  specimen.  Hence  the  specimen 
is  loaded  only  once.  A  strain  gage  mounted  at  the  middle  of  the  output  bar  is 
used  to  monitor  the  strain  pulses,  in  particular  the  first  transmitted  pulse,  in 
the  output  bar. 


5.2  Specimen  Preparation 


In  the  experiments  a  single  SMA  specimen  345  mm  long  was  used  as  well 
as  two  short  specimens  of  25.4  mm  length.  All  the  specimens  had  diameter 
12.7  mm.  After  machining  the  specimens  to  the  appropriate  lengths  they  were 
heated  to  540  in  standard  atmosphere  for  2  hours  and  furnace  cooled.  This 
process  was  used  to  erase  history  of  prior  plastic  deformation.  A  thin  oxide 
layer  was  formed  during  the  heat  treatment,  but  this  did  not  affect  the  overall 
response  of  the  material.  In  the  long  bar,  six  strain  gauges  were  placed  at  dis¬ 
tances  10  mm,  20  mm,  40  mm,  80  mm,  160  mm  and  320  mm  from  the  impact 
end.  A  high  temperature  strain  gauge  adhesive  was  used  and  the  specimens 
were  then  annealed  at  100  C  for  1  hour.  Subsequently,  the  specimens  were 
cooled  to  dry  ice  temperature  (-70  °C)  and  then  brought  to  room  temper¬ 
ature  for  testing.  All  tests  were  performed  at  room  temperature  (nominally 
20  °C).  A  Differential  Scanning  Calorimeter  (DSC)  was  used  to  determine  the 
transformation  temperatures  in  the  material.  As  can  be  seen  from  the  DSC 
measurements  shown  in  Figure  14,  under  the  indicated  temperature  cycling, 
the  specimens  were  in  a  twinned  martensitic  state  during  the  tests.  In  order 
to  obtain  preliminary  information  on  the  mechanical  behavior  of  this  material 
quasi-static  compression  test  was  performed  on  one  of  the  short  specimen  in  a 
standard  testing  machine.  Since  the  dynamic  test  involved  only  detwinning  of 
martensite  the  quasi-static  tests  were  done  at  room  temperature.  These  tests 
were  used  to  obtain  the  stiffness  of  the  martensitic  phase  and  the  critical 
stresses  and  u/  for  onset  and  finish  of  detwinning. 

The  material  constants  used  for  the  detwinning  model  are  summarized  in 
Table  3.  The  hysteresis  simulated  by  the  model  (Section  2.3)  and  the  actual 
hysteresis  from  the  quasi-static  test  are  given  in  Figure  15. 
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Table  3 


Material  parameters  used  in  the  SMA  model  for  det winning 


Material  constant 

Value 

Description 

EM 

42  X  10^  Pa 

Modulus  of  elasticity  in  martensite 

H 

0.027 

Maximum  detwinning  strain 

0-5 

-125MPa 

Start  of  -4  M'^  deformation 

-273MPa 

Completion  of  M*  deformation 

Y^d 

CsH 

pb^ 

Tempersture  (’C) 

Fig.  14.  Differential  Scanning  Calorimeter  measurements  of  the  SMA  specimen. 
5.3  Dynamic  Results 


As  indicated  earlier,  in  the  Hopkinson  bar  experiment  a  345  mm  long  rod  in¬ 
strumented  with  six  strain  gauges  was  placed  behind  the  input  bar.  The  output 
from  these  gauges  is  shown  in  Figure  16.  Strain  gauge  number  3  (40mm)  suf¬ 
fered  a  partial  debond  during  the  test  and  hence  the  results  from  this  gauge 
are  not  meaningful  beyond  the  point  marked  by  the  dark  dot  in  the  figure. 
The  elastic  wave  in  the  input  bar  was  not  recorded  due  to  an  error  in  the  de¬ 
vice;  all  other  gauges  worked  well  and  recorded  the  strain  profile  as  the  wave 
propagated  down  the  length  of  the  SM.4  rod.  An  x-t  diagram  corresponding 
to  elastic  wave  propagation  in  this  specimen  is  shown  in  Figure  17.  The  strain 
gauge  locations  are  indicated  by  the  thin  vertical  lines  and  the  leading  edge 
of  the  initial  loading  pulse  is  shown  by  the  dark  line;  this  pulse  reaches  each 
one  of  the  gauges  at  the  time  where  the  dark  line  intersects  the  vertical  lines. 
From  the  timing  of  the  elastic  wave  arrival  at  each  gauge,  the  elastic  wave 
speed  was  determined  to  be  2500  m/s.  The  elastic  wave  reaches  the  far  end 
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Fig.  15.  Quasi-static  hysteresis  of  SMA  specimen  and  the  model  simulation. 


Fig.  16.  Strains  measured  by  the  gauges  mounted  on  the  SMA  bar.  Gauge  3  suffered 
a  partial  debond  at  the  point  indicated  by  the  dark  circle  and  hence  the  data  beyond 
this  time  should  not  be  interpreted. 

of  the  specimen  about  138  /rs  after  impact.  The  duration  of  the  loading  pulse 
is  about  90  ps  and  hence  an  unloading  pulse  propagates  from  down  the  spec¬ 
imen  with  the  elastic  wave  speed  (since  the  unloading  is  elastic).  This  wave 
is  shown  by  the  line  with  an  arrow  at  the  tip.  Time  i  =  0  corresponds  to  the 
first  arrival  of  the  loading  pulse  at  the  strain  gauge  in  the  input  bar. 

As  seen  in  Figure  16,  the  strain  in  the  first  two  gauges  increases  rapidly  to 
a  level  of  about  1.3%  and  levels  off  as  the  load  from  the  input  bar  levelled 
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Fig.  17.  X-t  diagram  indicating  arrival  of  the  elastic  wave  front  at  the  gauge  loca¬ 
tions. 

off.  The  oscillations  seen  in  these  gauges  near  the  plateau  are  Pochhammer- 
Chree  oscillations  that  appear  in  bars.  At  around  290  pLS  the  unloading  wave 
from  the  end  of  the  loading  pulse  reaches  the  first  two  gauges  and  the  strain 
begins  to  decrease;  however,  because  the  strains  beyond  0.3%  were  the  result 
of  detwinning  (see  the  quasi-static  results  in  Figure  15),  these  strains  are  not 
recovered  and  a  permanent  strain  of  about  1%  is  left  at  these  locations.  The 
signal  in  gauge  4  clearly  indicates  the  dispersion  of  the  wave  -  higher  strain 
levels  propagate  at  significantly  slower  speeds  and  arrive  later  at  the  gauge 
location.  Hence  a  broadening  of  the  strain  pulse  can  be  seen  -  the  peak  in 
the  strain  at  gauge  4  occurs  7b  i^s  after  elastic  wave  arrival  while  it  occurs  in 
about  20  fis  in  gauge  1.  This  delay  also  results  in  the  peak  strain  not  being- 
sustained  for  too  long  as  the  elastic  unloading  pulse  reaches  the  gauge  quickly; 
once  again  a  residual  strain  of  peak  strain  -  0.3%  is  left  at  this  gauge  location. 
The  same  behavior  is  seen  in  gauge  5  where  due  to  its  distance  from  the  impact 
end,  and  due  to  the  slowness  of  the  inelastic  waves,  the  peak  strain  reached  is 
only  about  0.5%.  Once  again  a  residual  strain  is  left  in  this  location.  In  gauge 
6,  the  reflected  wave  from  the  end  of  the  SMA  rod  (left  free  in  this  experiment) 
causes  unloading  of  the  gauge;  a  very  small,  but  measurable  permanent  strain 
or  detwinning  is  observed  in  this  location.  Subsequent  to  the  test,  the  rod  was 
heat  treated  through  a  temperature  cycle  taking  it  above  A°^  first,  holding  for 
1  hour  and  then  cooling  below  M°7  and  warming  back  to  room  temperature. 
All  strain  gauges  recovered  their  original  state  indicating  full  recovery  of  the 
specimen. 
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The  results  of  this  experiment  can  be  used  to  extract  the  dynamic  stress- 
strain  response  by  applying  the  theory  of  one-dimensional  wave  propagation 
in  plastic  rods  due  to  (Rakhmatulin,  1945;  Von  Karman  and  Duwez,  1950; 
Taylor,  1958).  The  idea  is  a  simple  extension  of  the  rod  theory  for  elastic 
waves.  Let  us  assume  that  stress  is  only  a  function  of  strain,  i.e.  a  =  a(g:). 
Then  the  balance  of  linear  momentum  (1)  can  be  written  in  the  form 


Wft 


(55) 


Note  that  this  is  not  an  incremental  theory,  but  a  total  strain  theory;  therefore 
unloading  cannot  be  considered  here.  The  wave  speed  C{e)  of  disturbances  is 
no  longer  a  constant  as  in  the  linear  elastic  case,  but  a  function  of  strain: 


C{e)  =  (66) 

The  main  result  of  this  one  dimensional  theory  is  that  a  given  strain  (or  stress) 
level  will  propagate  into  the  rod  with  a  characteristic  speed  given  by  equation 
(56).  If  the  propagation  speed  of  strain  waves  in  a  one-dimensional  rod  is 
known  (measured  with  strain  gauges  as  in  the  experiment  discussed  above), 
equation  (56)  can  be  inverted  to  determine  the  stress-strain  behavior  of  the 
material:  ^  ^ 

<^(s)  =  /  =  p  /  C\Q<K  (57) 

This  representation  of  the  wave  speed  is  used  to  extract  the  constitutive  be¬ 
havior  of  the  material  (Bell,  1960;  Kolsky  and  Douch,  1962).  There  exists  a 
critical  point  in  the  stress-strain  curve:  o'{e)  —  0.  Strain  amplitudes  larger 
than  this  cannot  propagate  through  the  material.  Of  course,  in  the  experi¬ 
ment  discussed  above,  we  have  not  reached  this  stage;  in  fact,  this  would  be 
of  interest  in  determining  the  propagation  of  phase  transformation  fronts  and 
such  experiments  are  in  progress. 

The  propagation  speeds  of  different  strain  levels  were  obtained  from  the  results 
shown  in  Figure  16.  The  time  of  arrival  of  different  strain  levels  at  each  one  of 
the  five  gauges  were  determined  from  the  strain  measurements.  The  speed  of 
each  strain  level  C{e)  was  then  determined  from  the  known  distances  between 
the  gauges.  The  variation  of  the  wave  speed  with  strain  level  is  shown  in 
Figure  18;  a  smooth  trendline  is  also  shown  in  the  figure.  The  elastic  wave 
speed  is  about  2500  m/s  and  all  strain  levels  below  about  0.1%  travel  with 
this  speed;  this  suggests  that  there  is  really  no  significant  elastic  region  and 
that  even  small  strain  levels  are  susceptible  to  dispersion.  A  large  change 
in  the  wave  speed  occurs  at  around  0.3%  strain  which  corresponds  to  the 
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Fig.  18.  Variation  of  the  wave  speed  with  the  strain  level  as  determined  from  the 
strain  measurements.  The  line  is  an  eyeball  fit  to  indicate  the  data  trend.  Cubic  fits 
over  short  segments  were  used  to  determine  the  wave  speed  corresponding  to  each 
strain  level  in  the  determination  of  the  stress-strain  behavior. 


Fig.  19.  Stress-strain  response  evaluated  from  one-dimensional  wave  propagation 
measurements  (scatter  plot)  and  quasi-static  data  (solid  line). 

onset  of  massive  detwinning  deformation.  Beyond  this  level,  the  wave  speed 
drops  to  about  1000  m/s  and  varies  more  slowly.  If  the  averaged  data  on  the 
wave  speed  variation  with  strain  level  is  used  in  equation  (57),  the  resulting 
numerical  integration  provides  the  stress  strain  relationship  associated  with 
the  detwinning  deformation  in  the  SMA  rod.  Such  a  relationship  is  shown  in 
Figure  19.  The  scatter  in  the  plot  is  a  result  of  the  averaging  of  the  noisy  data 
in  Figure  18;  the  solid  line  shows  the  trend  of  the  data. 
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5-4  Numerical  simulations  of  the  Hopkinson  bar  experiment 


A  numerical  simulation  was  performed  and  results  were  compared  with  the 
experimental  data.  As  indicated  earlier,  due  to  a  trigger  failure,  the  signal 
in  the  input  bar  was  lost  so  only  the  readings  of  the  six  strain  gauges  on 
the  specimen  were  available.  In  order  to  supply  proper  boundary  conditions 
the  signal  from  the  first  strain  gauge  (at  10mm)  was  used  and  the  remaining 
gauges  were  simulated.  Gauge  number  3  was  not  included  in  the  modelling 
because  it  unglued  during  the  test. 

The  Hopkinson  bar  experiment  was  done  at  room  temperature  and  due  to 
the  heat  treatment  of  the  specimen  prior  to  the  test  it  was  in  fully  twinned 
martensitic  state.  The  SMA  model  was  applied  in  detwinning  conditions  (Sec¬ 
tion  2.3)  with  the  material  constants  given  in  Table  3.  The  measured  stress- 
strain  response  at  room  temperature  and  the  simulated  hysteresis  are  shown  in 
Figure  15.  The  adaptive  FEM  scheme  was  chosen  because  of  its  accuracy  and 
ability  to  predict  precisely  the  positions  of  the  both  elastic  and  transformation 
shocks.  The  results  are  presented  in  Figure  20. 

As  expected  from  the  numerical  examples  studied  in  section  4.1  the  strain 
wave  splits  into  an  elastic  and  a  transformation  front.  The  transformation 
front  timing  and  magnitude  at  all  strain  gauges  is  in  good  agreement  with  the 
experiment.  The  small  oscillations  observed  at  the  first  two  gauges  are  due  to 
surface  effects  caused  by  the  impacting  projectile.  Such  effects  cannot  possibly 
be  modelled  within  a  1-D  formulation. 

There  is,  however,  a  noticeable  disagreement  in  the  timing  of  the  elastic  fronts. 
The  reason  for  this  is  the  deviations  from  linear  behavior  for  small  strains.  The 
polynomial  model  always  predicts  a  linear  response  until  the  beginning  of  the 
detwinning  deformation.  However  an  inspection  of  Figure  15  shows  a  smoother 
nonlinear  stress-strain  relationship  for  small  strain  values. 

To  verify  the  hypothesis  that  the  disagreement  is  due  to  the  initial  elastic 
response  of  the  model  an  independent  numerical  simulation  of  the  dynamic 
experiment  was  performed.  A  phenomenological  deformation  plasticity  model 
was  used  instead  of  the  constitutive  model  of  Section  2.3.  The  loading  is  as¬ 
sumed  to  have  the  form  of  a  sixth  degree  polynomial  that  curve  fits  the  loading 
part  of  the  quasi-static  hysteresis  in  Figure  15.  The  unloading  was  assumed 
linear,  the  slope  being  the  modulus  of  martensite,  A2GPa,  as  measured  by 
the  quasi-static  experiments.  Due  to  the  fact  that  the  deformation  is  mostly 
detwinning  of  martensite  there  is  no  significant  release  of  latent  heat,  so  the 
quasi-static  hysteresis  is  very  close  to  the  actual  material  behavior  in  the  dy¬ 
namic  case  (Figure  19). 

The  results  of  the  simulation  of  the  dynamic  problem  are  shown  in  Figure  21. 
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Fig.  20.  And  adaptive  FE  analysis  of  experimental  data  under  isothermal  conditions. 
The  first  strain  gauge  is  used  to  define  the  boundary  condition. 


Time  (microseconds) 


Fig.  21.  And  adaptive  FE  analysis  of  experimental  data  under  isothermal  conditions 
and  a  curvefit  of  the  hysteresis.  The  first  strain  gauge  is  used  to  define  the  boundary 
condition. 

This  time  the  wave  profiles  are  matched  much  more  closely  and  the  small 
disagreements  can  be  attributed  to  measurement  errors  and  effects  of  lateral 
inertia  not  included  in  the  simulation.  It  should  be  noted  that  unlike  a  consti¬ 
tutive  model  based  on  physical  principles  such  an  approach  will  only  work  for 
a  particular  SMA  specimen  and  particular  operating  temperature.  However 
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using  a  curve  fit  for  the  loading  part  of  the  hysteresis  is  sufficient  to  check 
whether  disagreements  between  experiments  and  simulations  are  indeed  due 
to  the  constitutive  model. 


6  Conclusions 


The  problem  of  dynamic  loading  of  one-dimensional  polycrystalline  SMA  rods 
has  been  explored  numerically  and  experimentally.  FEM  simulations  were  per¬ 
formed  for  SMAs  experiencing  both  pseudoelastic  phase  transformation  as  well 
as  detwinning  deformations.  A  long  SMA  rod  was  tested  in  a  split  Hopkin- 
son  bar  experiment  under  detwinning  conditions.  The  wave  propagation  was 
observed  through  strain  gauges  placed  on  the  surface  of  the  specimen.  The 
strain  history  at  the  gauge  locations  obtained  through  numerical  simulations 
of  the  dynamic  experiment  was  found  to  be  in  good  agreement  with  the  actual 
results. 

Computational  solutions  were  shown  to  coincide  with  known  analytical  re¬ 
sults.  Nonlinear  shock  formation  and  velocities  were  captured  correctly  by  the 
FEM  simulations.  The  standard  semi-discrete  FEM  approach  for  hyperbolic 
problems  was  complemented  by  an  adaptive  mesh  refinement  technique.  The 
utilization  of  the  Zienkiewicz-Zhu  error  indicator  lead  to  an  order  of  mag¬ 
nitude  decrease  of  the  computational  time.  Energy  dissipation  calculations 
for  both  detwinning  of  martensite  and  stress-induced  phase  transformation 
showed  that  the  strain  energy  can  be  reduced  by  80-90%  which  suggests  that 
SMAs  can  be  used  effectively  as  shock-absorption  devices. 

On  the  experimental  side,  it  has  been  shown  that  an  instrumented  Hopkinson 
bar  can  be  used  effectively  to  evaluate  the  wave  and  phase  propagation  char¬ 
acteristics  in  the  SMA  rods.  Through  the  use  of  multiple  strain  gauges,  the 
phase  velocity  at  the  different  strain  levels  was  obtained  easily.  An  inelastic 
deformation  theory  was  used  to  interpret  the  dispersion  in  terms  of  the  un¬ 
derlying  dynamic  material  response  of  the  material.  Dynamic  and  quasi-static 
material  response  were  shown  to  be  in  excellent  agreement. 

Through  careful  calibration  of  the  constitutive  model  for  SMAs  the  peak  strain 
levels  of  the  Hopkinson  bar  experiment  were  accurately  predicted.  The  main 
drawback  of  this  model  is  its  initial  linear  response  in  the  case  of  detwinning 
and  the  existence  of  kinks  in  the  hysteresis  curve.  Accurate  predictions  of  the 
entire  experimental  data  were  obtained  by  using  a  polynomial  curve  fit  of  the 
quasi-static  hysteresis  of  the  material.  Both  the  wave  timings,  shape  and  peaks 
were  modelled  within  experimental  error. 

The  material  and  environmental  conditions  used  in  the  Hopkinson  bar  exper- 
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iments  correspond  to  a  detwinning  deformation  of  the  martensitic  phase,  but 
the  methods  can  be  easily  adapted  to  stress  induced  martensitic  transforma¬ 
tion  in  tests  at  higher  temperatures.  Theoretical  work  can  also  be  extended 
to  more  realistic  2-D  and  3-D  geometries.  Complicated  SMA  components  and 
structures  can  be  simulated  to  better  understand  the  nonlinear  wave  propa¬ 
gation  phenomena  as  well  as  the  practical  aspects  of  their  energy  dissipation 
capabilities.  More  refined  models  which  incorporate  both  detwinning  and  pseu¬ 
doelastic  deformations  simultaneously  and  also  predict  accurately  the  smooth 
hysteresis  of  the  detwinning  deformation  will  be  extremely  helpful  in  further 
studies  of  wave  propagations  in  polycrystalline  SMAs  and  are  currently  under 
consideration. 
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